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A leading approach to the modeUing of extreme mass ratio inspirals involves the treatment of 
the smaller mass as a point particle and the computation of a regularized self-force acting on that 
particle. In turn, this computation requires knowledge of the regularized retarded field generated 
by the particle. A direct calculation of this regularized field may be achieved by replacing the point 
particle with an effective source and solving directly a wave equation for the regularized field. This 
has the advantage that all quantities are finite and require no further regularization. In this work, 
we present a method for computing an effective source which is finite and continuous everywhere, 
and which is valid for a scalar point particle in arbitrary geodesic motion in an arbitrary background 
spacetime. We explain in detail various technical and practical considerations that underlie its use 
in several numerical self-force calculations. We consider as examples the cases of a particle in a 
circular orbit about Schwarzschild and Kerr black holes, and also the case of a particle following a 
generic time-like geodesic about a highly spinning Kerr black hole. We provide numerical C code 
for computing an effective source for various orbital configurations about Schwarzschild and Kerr 
black holes. 



I. INTRODUCTION 

There has been much recent interest in the study of Extreme Mass Ratio Inspiral (EMRI) systems. These systems 
typically involve a compact, solar mass object inspiralling into an approximately million solar mass black hole. Such 
massive black holes are expected to exist at the center of most galaxies Ij. 

EMRIs are expected to provide a strong source of gravitational waves for future generations of gravitational wave 
detectors [2H1]. There is also hope that parameters for these sources can be accurately estimated, enabling studies 
and measurements of the strong field region of central supermassive black holes [SH?]. In order to achieve accurate 
parameter estimation, it is essential that highly accurate gravitational waveforms are available. This, in turn requires 
highly accurate, long-time models of the inspiral. 

A leading approach to the accurate modelling of EMRI systems arises from the fact that the mass ratio, fi, is very 
small. This makes it possible to treat the system within perturbation theory, in which the smaller object is assumed 
to be a point particle generating a perturbation about the background of the larger mass. At zeroth order in /x, the 
smaller object merely follows a geodesic of the background. At first order, it deviates from this geodesic due to its 
interaction with its self-field. This deviation may be viewed as a force acting on the smaller object, referred to as the 
self-force. The calculation of this self-force is critical to the accurate modelling of the evolution of the system. 

A naive calculation of the first order perturbation leads to a retarded field which diverges at the location of the 
particle. The self-force, being the derivative of the field, therefore also diverges at the location of the particle and 
must be regularized. A series of derivations of the regularized first order equations of motion (now commonly referred 
to as the MiSaTaQuWa equations, named after Mino, Sasaki, Tanaka [S] and Quinn and Wald [1] who first derived 
them) for a point particle in curved spacetime have been developed [8lll7j. culminating in a recent rigorous work by 
Gralla and Wald [TH] and Pound [19J in the gravitational case and by Gralla, et al. [2D] in the electromagnetic case. 
Several practical computational strategies have developed from these formal derivations: 

• By measuring the flux of gravitational waves onto the horizon of the larger black hole and out to infinity, a 
time-averaged dissipative component of the self-force - which is finite and does not require regularization - may 
be computed. This, however, neglects potentially important conservative effects which may significantly alter 
the orbital phase of the system. 

• The mode-sum approach, introduced in Refs. |2 H 122 ) . which involves the decomposition of the retarded field 
into spherical harmonic modes (which are finite, but not differentiable at the particle), solving for each mode 
independently and subtracting "regularization parameters" , then summing over modes. This method has been 
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used to compute the self force for a variety of configurations in the Schwarzschild and Kerr [551411] 

spacetimes. 

• The effective source approach |42fl47j in which the regularization is done before solving the wave equation. In 
this case, all quantities are finite throughout the calculation and one directly solves a wave equation for the 
regularized field. A review of this approach can be found in Note that the effective source proposed by 
Lousto and Nakano 03 differs in that it is not derived from the Detweiler- Whiting singular field. 

• The matched expansion approach j49l I50j in which a quasi-local expansion of the Green function |51fl53j (which 
is valid in the recent past) is matched onto a quasi-normal mode sum (valid in the distant past|^ The retarded 
field is then computed as the integral of this matched retarded Green function along the worldline of the particle. 

For a comprehensive review of the self-force problem, see Refs. [551 - f57] . The present work focuses on the third of these 
strategies, the effective source approach. In this approach, the point particle source is replaced with a finite effective 
source leading to a wave equation which admits the correct regularized field (at the particle) as a solution. 

Given our motivation in studying the EMRI problem, it is the gravitational self- force which is of the most interest. In 
this paper, however, we instead study the analogous scalar self-force. This allows us to develop insight and techniques 
without being obscured by the additional complexity of the gravitational case. It should be noted, however, that this 
extra complexity is predominantly only calculational and comes in the form of larger expressions. Conceptually, the 
calculations done here follow through for the gravitational case with few modifications. 

The purpose of this paper is to provide a comprehensive exposition of the scalar effective source employed in a 
variety of recent and ongoing self- force calculations 43, 47l[58l[59]. Starting with its covariant definition, we present 
its coordinate construction and the various modifications that we found necessary in order to get the effective source 
to its current "best" form. Much of the paper is technical in nature, but we believe that all the details provided here 
are essential to anyone interested in pursuing an effective source approach to self-force calculations. 

The layout of the paper is as follows. In Sec. [lljwe introduce the effective source approach in detail and compute 
approximations to the singular field and effective source in the form of covariant expansions. In Sec. |III[ we develop 
practical methods for evaluating these approximations in a specific spacetime in terms of coordinate expansions. We 
give example calculations for the case of a circular geodesic orbit in Schwarzschild and Kerr spacetimes and a generic 
orbit in Kerr spacetime in Sec. |IV| In Sec. |V] we conclude with a discussion on aspects of the calculation and on 
prospects for future applications. In Appendix |Xj we develop covariant expansions of various biscalars used in this 
paper. In Appendix [B| w e discuss a modification to the covariant expansion which yields substantial practical benefits. 
Finally, in Appendix |C[ we discuss issues related to efficient numerical implementations for computing the singular 
field and effective source. 

Many of the expressions developed in this work, although useful, are too unwieldy to be given in printed form. 
Instead, we have made available all expressions we deem to be useful online [30] as Mathematica code. Furthermore, as 
this work is intended to provide computational tools for those interested in doing self- force calculations, we also include 
a library of C code for computing the singular field and effective source for various configurations in Schwarzschild 
and Kerr spacetimes. The intention is for this code to be a "black box" which can be easily incorporated into existing 
numerical codes, whether they are 3 -I- ID, 2 -t- ID or 1 + ID. 

Throughout this paper, we use units in which G — c — 1 and adopt the sign conventions of |61j . We denote 
symmetrization of indices using brackets (e.g. (a/3)) and exclude indices from symmetrization by surrounding them 
by vertical bars (e.g. (a|/3|7)). Roman letters are used for free indices and Greek letters for indices summed over 
all spacetime dimensions. Roman letters starting from i are used for indices summed only over spatial dimensions. 
Capital letters are used to denote the spinorial/tensorial indices appropriate to the field being considered. For 
convenience, we frequently make use of the shorthand notation of Ref. [27] by introducing definitions such as Rucma\a = 



II. EFFECTIVE SOURCE APPROACH 



To compute the self- force, acting on a point particle with scalar charge g, knowledge of the retarded field, 'fret, 
generated by the particle is required. This field is a solution of the inhomogeneous wave equation, 

5^{x 



z{t)) 



cLt, 



-g 



(1) 



^ In black hole spacetimes there is also a branch cut integral which must be evaluated in the region where the quasi-normal mode sum is 
used. Substantial recent progress has been made towards the calculation of this branch cut contribution I54| . 
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where the source corresponds to a point particle on a worldhne 7 in some background spacetime and where 

2? = (□ - ^R) (2) 

is the scalar wave operator. A naive calculation of the self-force from this retarded field will diverge when evaluated 
at the location of the particle. In order to compute a meaningful self-force, one must therefore find a regularized 
retarded field. This may be achieved by separating the field into singular (S) and regular (R) parts, 

«'rct = *s + $R- (3) 

The identification of a singular field which gives the correct regularized self-force is crucial. Using a Green function 
decomposition, Detweiler and Whiting |14j were able to find a representation of the singular field which is valid in a 
region near to the particle. It is a solution of the same inhomogeneous wave equation ([l]) as the retarded field. A 
brief overview of their approach is given in the next subsection. 

Given knowledge of the singular field, one must then prescribe a method of computing the regularized field. In the 
effective source approach, first proposed independently by Barack and Golbourn [33] and by Vega and Detweiler [35] , 
the splitting of the self-field into regular and singular parts is done at the level of the wave equation, 

2?$ret, - + 2?*R, (4) 

before solving for the field. One then solves directly the equation for the regularized field. 

I?$R = -47r<7 / ^!(^^iMld7- _ ^ 0. (5) 



The regularized self-force is then simply given by the derivative of this regularized field, 

r = gv-^R. (6) 

This method has several advantages: 

• It does not rely on the separability of the field equations. This is particularly important in the Kerr spacetime 
where the perturbation equations are not fully separable in the time domain. 

• There are no troublesome delta functions or singularities to deal with. This is particularly advantageous in 
numerical calculations where smoothness is desirable. 

• In comparison to methods which first compute ^rct and then regularize, there is no need to cancel two large 
quantities ($3 and $rot) to get the self- force, so in principle the field one solves for is inherently more accurate. 
This is particularly relevant in numerical calculations, where the cancellation of large quantities may lead to 
considerable round-off errors. 

• Its applicability in the time domain means that the orbit may be evolved, coupling the geodesic equations into 
the wave equation and source calculation. 

In principle, the regularized field is a solution of the homogeneous wave equation and the self-force is determined 
purely from the boundary condition for <I>r = <i>i.ct — ^s- However, in practice the singular field identified by Detweiler 
and Whiting is not defined globally (it is not even clear that a global definition exists) . One must therefore introduce 
a method for restricting the singular field to a region near the particle. Furthermore, in practice an exact calculation 
of the singular field away from the particle proves difficult; it is much easier to calculate an approximation to the 
singular field, denoted by $s and to solve for an approximate regularized field ^r,. The construction of an approximate 
singular field must ensure that its local expansion near the particle matches that of the actual singular field sufficiently 
well that evaluating the self-force using $r yields the correct value at the particle. It is important to note, however, 
that this approximate regularized field becomes meaningless far from the particle. 

There are two different approaches to dealing with the problem of the lack of a global definition for the singular 
field. Vega and Detweiler [35] tackle the issue of restricting the singular field to a region near the particle with the 
use of a window function, W, and split the retarded field as 

$rot = W^s + $R- (7) 

The window function is chosen so that near to the particle W^s remains a good approximation to the singular field, 
while far away from the particle W dies away sufficiently quickly that $r w $rct- Barack and Golbourn [33] take an 
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alternative approach. They introduce a world-tube around the particle. Inside the world-tube, they solve for $p{ and 
outside they solve for $rot- They then impose ([3| as what is essentially a "change of variables" on the world-tube 
boundarjj^ In this case, the lack of a global definition for $g is no longer an issue as the only requirement on l>g is 
that it approximates the singular field sufficiently well near the particle. 

In both cases, the approximation to the singular field, $Sj is no longer a solution of Eq. ([T]). The source term now 
has additional structure away from the particle, extending throughout the worldtube or the region of support of the 
window function. As a result, the approximate regularized field is now a solution of the inhomogeneous wave equation 
with an effective source, S'cff : 

V^n = -4^9 / ^-li^-^^dr - V (w^s) = S.g. (8) 
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In contrast to Eq. ([ij, however, this source has the advantage of being regular and smooth everywhere except at 
the location of the particle where it is still regular, but of finite differentiability, the level of differentiability being 
determined by the choice of approximation to the singular field. 



A. Exact expression for the singular field 

In order to obtain an expression for the singular field, we follow Detweiler and Whiting [W in introducing the 
Hadamard form ^62. ,63] for the singular Green function, 

Gs(x, x') = i [C/(x, x')d{a{x, x')) + V{x, x')0{a{x, x'))] , (9) 

which is obtained by adding a homogeneous solution (in this case V{x,x')) of the wave equation to the symmetric 
Green function, Gsym = 5 (Grct + Gadv) [M]- Here, 6{a{x,x')) is the covariant form of the Dirac delta function, 
9 {a {x, x')) is the Heaviside step function, and U{x,x') and V{x,x') are symmetric biscalars which are regular for 
x' — X. The biscalar a {x, x') is the Synge |551 165j world function, which is equal to one half of the squared geodesic 
distance between x and x' . 

This singular Green function is a solution of the same wave equation as the symmetric Green function, but differs 
in that it has support only on and outside the light-cone. Note that this singular Green function is not guaranteed to 
exist globally. Its definition depends on the existence of the unique function V{x,x'), which is only true provided x 
and x' are within a convex normal neighborhoocj^ Fortunately, in the effective source approach we only require that 
it exists in a neighborhood of the particle, in which case it can be given the clear definition We now define the 
singular field by 

$s= / Gs{x,z{T))dT. (10) 



Substituting ^ into ( 10 1 and making the change of variables r — >■ a{x, z{t)), we obtain an expression for the singular 
field which depends on a finite portion of the particle's world line: 

^ , , U{x,x') U{x,x") 1 f\^, , 

<^g(x) ^ — f _ — Ir + - j V{x,z{T))dT. (11) 

Here we have introduced the retarded and advanced points x' and x" corresponding to the retarded and advanced 
times u and v on the world-line 7 associated with the field point x (Fig.[T]). 



B. Approximation to the singular field 



The expression for the singular field given in Eq. ( 11 ) is very general. It is valid for any worldline in any spacetime 
provided the field point x is sufficiently close to the worldline that the singular Green function can be defined. In 



^ In a numerical implementation, it is common to reduce the wave equation to a system of first order equations. In this case, it may be 

necessary to impose the conditions not only on <I>i^, but also on its derivatives. 
^ When considering Hadamard form Green functions such as |9]l, one typically defines them within a causal domain |63j . The singular 

Green function is acausal so this must be relaxed to a definition within a convex normal neighborhood, requiring only that (t(x, x') be 

unique. 
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FIG. 1. The singular field at the point x can be expressed in terms of the retarded and advanced distances to the world-line, 7. 



practice, it is only in very simple spacetimes that U{x, z{t)), cr(a;,z(r)) and V{x,z{t)) may be computed exactly. 
In many curved spacetimes of interest (including Schwarzschild and Kerr) this is not the case and one must find an 
approximation to (11). 



In the present work, we choose a covariant series expansion of ( 11 ) (taken to second order in the geodesic distance 
from the field point to the world line) as a starting point for our approximation to $s- In doing so, we use the methods 
described in Refs. [SS] and [57] to consolidate the dependence of $3 on the advanced and retarded points x' and x" 
into a single arbitrary point x on the worldline. This has the additional advantage of making the dependence of x'{x) 
and x"{x) on x explicit, so that x is truly an arbitrary point on the worldline (sufficiently close to x' and x") with no 
implicit dependence on x. We additionally make use of the techniques of Ref. [52j to compute covariant expansions 
of all required bitensors. 

Given the primary motivation of studying black hole spacetimes such as Schwarzschild and Kerr, it is reasonable to 
assume that the spacetime is vacuum (i.e. Rab = 0). However, in the present work, we do not make that assumption. 
This is motivated by the fact that some of the leading-order terms in the covariant local expansion of the gravitational 
singular field involve the Riemann tensor and do not vanish in vacuum. In contrast, the analogous leading-order terms 
for the scalar case involve only the Ricci tensoij^ As far as the covariant local expansions are concerned then, the 
scalar singular field in a non-vacuum spacetime best captures the structure of the gravitational singular field in a 
generic spacetime. By not assuming vacuum, we are therefore emulating some of the extra complexity which would 
otherwise only appear in the gravitational case. 

The covariant expansion of <I>s requires, in turn, the expansion of the functions U{x,x'), U{x,x"), Ua'u'^ , (Tq"u" 
and V{x, z{t)) about the point x. We compute these expansions in Appendix|Aj Substituting ( A9), ( AlO), (All ) and 
(|A16|) into (fnl), we get 



6s3 



12 s 



2ri?„ 



Rn 



Ruu ( ^ 



^--^}R-s 



24s 



Racr\cr + {Racr\u ~ '^Ru(j\a) f + (2i?«o-|u - Ruu\a){y'^ + ^) + Ruu\u{y'^ + 3 rS^)) 



+ 4(e-g)(i?|«rs-i?.s)+2^i 



(12) 



o-gM" (the projection 



where s = (5"" -f u^u' )a-aCFfj (i.e. the projection of aa orthogonal to the worldline), and r : 
along the worldline) and we adopt the notation of Haas and Poisson [27] in defining -R„ct«o- la- 
Letting e be a measure of the geodesic distance from x to the world- line (i.e. x), the first term here is 0{e~^), the 
second group of terms is 0{e^) and the third group of terms is 0{e'^). The difference between <I>s and $3 is then ©(e^) 



R^^^s.^.u^a^u^a'a'. 



* More specifically, in the scalar case the first four orders in the covariant expansion of the tail term, V{x,x'), involve only the Ricci 
tensor and Ricci scalar. As a result, V{x, x') = ©(e**) in vacuum and the tail term could be neglected in the present calculation. In the 
gravitational case, however, the tail term, V^^' bb' > ^' ) ^ h^-^ ^ leading order component involving the Riemann tensor. This means that 
even in vacuum Vi^^n,i,i (x, x') = 0(1). 
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C. Approximation to the effective source 

Given the approximation ( |12[ ) to the singular field, a corresponding effective source may be computed by applying 
the wave operator to i>s- (This requires cancelling the divergent terms in $5, so the derivatives in the wave operator 
must be computed very accurately. In particular, straightforward numerical differentiation does not provide sufficient 
accuracy close to the particle.) In this section, we give an exact expression for the effective source and compute 
an approximation which is valid near the particle. This approximation gives insight into the properties of a source 
derived from a particular order approximation to the singular field. 

Before proceeding further, we will clarify the meaning of 'order' as used in this context. All approximations are 
considered as expansions in powers of e, which is roughly speaking the distance between x and the world-line (i.e. the 
length of the bivector ct"). This means that s, r, and ct" are all of order e. The order of an approximation is then 
defined in terms of the order of the approximation to the singular field. The first order approximation is given by the 



leading term (of order e""'^) in the approximation to $g, i.e. the first term in (12). The second order approximation 
is given by the first two orders (to order e*') in the expansion of $s- As there is no term at order e" in ( [l2| ), at this 
stage the second order approximation is equivalent to the first order approximation. As will be discussed in Sec. |III| 
this will not, however, always be the case. Likewise, the third order approximation includes terms up to order e in 
$8 and the fourth order includes terms up to order e^. When referring to the effective source, the order referred to 
will be determined by the order of the singular field from which it is derived so that the first order effective source 
will be given by the wave operator acting on the first order singular field, and so on. 



1. First order 



The first order approximation to the singular field is given by the leading term in ( 12 ) 



(1) 



(13) 



which is of order e ^ . Since the wave operator contains second derivatives, one would in general expect that the result 
would be of order e~^. Applying ([2]) to ( 13 1 we find that this does appear to be the case: 



5(i)_p.l)_3(f^ 



— ^ VarV"r- — 



rDf 

e3 



(14) 



Noting that V qT and ctq." are 0(1) and □ r is 0(e), we see that the first three terms here are 0{e^^) and the last two 
are ©(e^^); it would appear that the first order effective source is 0{€^^). However, expanding Var, ctq", □ r and R 
about X, we find that the 0{e~^) terms cancel, leaving a source which is O(e-i): 



c.(l) _ 
''eft — 



3r2-s2 



3s5 
3r2 



R,. 



6 



R 



1 

it 

U(7U(7\(T 



3 (2ri?, 



Raa) 



12 S3 



[rR. 



(7a\u 



6rR 



'U(t\(T 



3R. 



(T(t\(J ) 



s 



0{e). 



(15) 



Here, the first group of terms are 0{e~^) and the second group are 0(1). In other words, the first order source 
diverges at the particle like 1/e, i.e. it is C^^ . This is sufficient to give a finite, but discontinuous regularized field. As 
a result of the discontinuity of the field at the particle, it is not possible to compute the self-force from its derivative. 

Note that the second order source will be the same as the first order source and will therefore have the same 
properties, with one caveat: by decomposing into m-modes, Barack et al. were able to extract a self-force from a 
second order source. This may be understood as a result of the "averaging" effect the Fourier transform used in the 
TO-mode decomposition has on the smoothness of the source. 



2. Third order 



The third order approximation to the singular field is given by the two leading terms in ( 12 ) 

I =2 



1 



-1 _q Ruaua ^ _ ( 2 r Rua 

OS'^ 12 



Raa + Ruu{r^ + S')) + ^ U - ^ I 



(16) 



7 



Applying ([2| gives the third order effective source: 



0(3) _ qW 

•^eff — ''off 



12 



rs 



(r □ r + Da) - (r^ + s^) - (f^ - s^) r V"f 



2i?(TaS 



-i + 6,e 



rs^Dr + s^DCT- (1 + V^r V"r)(3 - s 
+2i?„^s2 (^s2(s2 - r2)nr + r(3r2 -s^ -s2ncr + 3(r2 -s2))VarV"r 
+i?„„s2 r s2(3s2 - r2)n r + (f2 - s^) (3 + ~ s^Dcr + 3(7^ - s^) V^r V"f 
+ 2s4 Irs^a^'^p + 2{s^ - f2)a"^V'3r 

+2Ruaua 15?"* - 127^52 - 375^(72 - s2)n7 

-5^(372 - s2)ncr + 3(57* - 67^ + s'*)V„7V"7 



4R -5 



-2 -2 

r — s 



72- §2 



(17) 



which appears to have an additional contribution at 0{e^^) compared to the first order case. However, as before, 
re-expanding higher derivatives of a about x, we find that this exactly cancels the 0{e^^) contribution from the first 
order source, leaving a source which has a directional dependence at 0(1): 



0(3) _ c(l) 
''eff — ''off 



3s5 



3s= 



(2 rRua + Raa) — 



0{e) 



37^-5^ 
6s^ 



aua\a 



12s3 



(r R<T(j\u — 6 r Rua\a — 3-Rct<j|<j) 



0(e). 



(18) 



In this way, we see that including the third order contribution to the singular field gives a source which is C ^. This 
is now sufhcient to calculate both the regularized field and its derivative (i.e. the self- force). 



3. Fourth order 



Following the procedure once more, by including the fourth order contribution to $s, computing the associated 
effective source and re-expanding higher derivatives of a about x, we find that it has a contribution at 0(1) which 
exactly cancels that from the third order source, leaving a source which has a directional dependence at 0(e): 



:.(4) 



s. 



(3) 
off 



372 



6 



-R 



1 



uaua\a 



12 S3 



6rR 



'U(T\a 



3R 



a(T\a ) 



^R\c 



-0(e)=0(e). 



(19) 



Therefore, including the fourth order contribution to the singular field we obtain a source which is C'^. This is not 
only sufficient to give the self-force, but gives reasonably good convergence in numerical calculations. 



4- Higher orders 

We can clearly continue in this way (Fig. [2]), producing a smoother source at each step. Taking this expansion to 
its logical conclusion, if we can calculate $s exactly, then we find that 

= (20) 

and the self force comes purely from the boundary conditions. In practice, this would only require the computation 
of the expansion of $s to sufficiently high order to give an accurate numerical value in the region of interest. This 
may be difficult to impose in the window function approach, but in the world-tube approach the world-tube may be 
arbitrarily small and it may be possible. In that case, one would place a world-tube around the particle and then 
solve the system 



(21) 
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0{s-^) + 0{s) + Ois") 




□$sw 0(s-3) +0(s-i) + C(l) + 0{s) 



FIG. 2. Relation between order of approximation to the singular field and smoothness of the corresponding effective source. 
Also shown is the relation to the coefficients A^, B^, C^, D^, • ■ • in the large-Z expansion of the singular field used in the 
mode-sum regularization method [21] . 



where $i is the field inside the tube and $2 is the full retarded field outside the tube. The self force then comes from 
applying the change of variables (i.e., boundary condition) 

$2|w = $i|w + $slw (22) 

across the world tube boundary W. In this way, one may view the effective source as a correction for the fact that 
the singular field is not known exactly. 



III. COORDINATE EXPRESSIONS AND SOME PRACTICAL CONSIDERATIONS 



The previous section described the calculation of the singular field and effective source in a fully covariant manner. 
In practical applications, one needs to compute the singular field and effective source as a function of coordinate 
positions in a particular spacetime. A practical approach to doing so is to compute coordinate expansions of the 
singular field and corresponding effective source. In this section, we develop such expansions and give example 
implementations in Schwarzschild and Kerr spacetimes. In doing so, we exploit insight from the covariant approach 
to simplify the calculations as much as possible. 



A. Coordinate expansion of singular field 



All terms in Eq. (12 1 may be written in terms of erg and local quantities at x. In order to compute an explicit 
expression for a specific spacetime, it is convenient to expand Ua in the coordinate separation between x and x as 
follows [5T] : 



1. Write a{x,x) as a formal coordinate series expansion about x: 

1 



Ax"Ax^ + A^^^Ax'^Ax'^Ax^ 
+C^f,^s,Ax''Ax^Ax''Ax'^Ax' H 



B^^^sAx"Ax'^Ax^Ax^ 



(23) 



where Aa;° = — x° and where each of the coefficients is a function of x only and is symmetric in all indices. 

2. Differentiate this expression at x to get 

da =5a^Ax^ + (5,^^a + iA^^^)Ax^Ax^ + {A^^s^^ + AB^f^^s)Ax^ Ax^ Ax'^ 

+ (i?^^j,,a + 5C^^^s,)Ax^Ax^Ax'Ax' +■■■ (24) 

3. Use the identity 2a — (Jgcr" to recursively determine the coefficients Ag^^^, Bg^^^g. 



The result is a coordinate expansion of erg which may be substituted into (12). For the fourth order (i.e. C(e^)) 
approximation to the singular field, the coordinate expansion (23 1 must be computed to O \{Ax°')^]^ (i.e. the coefficients 
up to Cg^iyjj must be determined). Since the coefficients are just functions of the metric and its partial derivatives at 
a;, this calculation is easily achieved using a tensor software package such as GRTensorll ^67j or xCoba [551 IMj- Rather 
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than giving the full lengthy expressions, we present here only the leading two orders in the expansion in Schwarzschild 
spacetime to illustrate the structure: 

^ ' ' 2(f-2M) 2 2 ^ 2f 



. ^ ^ A A ^9 A A ,9 M 

2(f~2M)" 2 2 2f2 



2 Ar^ + -ArA6|2 + -ArA02 - — A^A^^ (25) 



where the point x is assumed to lie in the equatorial plane (the spherical symmetry of Schwarzschild means that it is 
always possible to ensure this is the case) . We provide a higher order expression for Kerr spacetime online [60^ . 
Next, we contract erg with the metric, Riemann tensor and four-velocity (all evaluated at x) to get r, s and Riemann 



terms such as Ruaua- We then substitute these into Eq. (12) to obtain the coordinate expansion of <i>s- In doing so. 



we only keep terms that contribute up to 0{e^). The first term in Eq. (12) is 0{e ^) and so requires the coordinate 



expansion of s to order O [(Aa;°)^] (equivalently, the expansion of to order O [(Aa;°)^] ). The second term is 0{e) 
and requires r^, and Ruaua to O [(Ax°)^] . The third term is ©(e^) and requires r^, and Ruaua\u to O [(Aa;'')^] 
(the leading order) and Ruaua\a to 0{0 [(Ax°)'^]) (again, the leading order). This results in an expression for the 
singular field which is valid to O [(Ax'')^] and has the general form 

~ ^ Q(2) + Q(3) + Q(4) + Q(5) , . 

' (&(2)+fc(3)+6(4)+&(5))3/2 ^ ^ 

where we use the notation a(„) = aa^...a^Ax°'^ ■ ■ ■ Ax"". 

Although there is a clearly defined 'true' singular field, in the effective source approach we may still view (&s as 
merely a computational tool with a certain degree of flexibility in choosing its particular form. Indeed, this coordinate 
approximation to the singular field is not unique - the only requirement it must satisfy is that it matches the 'true' 
singular field to a prescribed order - and it may therefore be replaced with any other expression which agrees with it 
to the desired order. 



The expression in the denominator of Eq. (26) is undesirable because it leads to long calculations, particularly 
when computing the derivatives required for the effective source corresponding to this choice of singular field. More 
importantly, the roots of this denominator are singularities in and, potentially, in S'eff. Since it is a power of 
a fifth-order polynomial, the denominator will have roots different from the trivial one. Ax" = 0, which represents 
the worldline of the particle. As a result, the effective source will have undesirable divergences at certain coordinate 
locations. (Note that it is C° at the location of the particle). Moreover, on any given time slice, the precise location 
of these singularities will depend sensitively on the position and four- velocity of the particle (on which the coefficients 
^Qia2 . ct„ depend). The presence of these extra singularities is purely an artifact of using a truncated series expansion 
to approximate s; the exact s increases monotonically away from the particle. This becomes problematic for any 
numerical application. 

For these reasons, it is advantageous to modify the singular field produced from the above described procedure. 
Noting that Ax"" = 0(e), we re-expand the coordinate expansion of $s about e = 0. In practice this is most easily 
achieved by introducing an explicit factor of e into the coordinate distances, Ax° eAx", expanding about e = 
(to 0{e^) for the fourth order singular field) and reading off the coefficient of each power of e. The result is an 
approximation to the singular field of the form 

^ C(6) + C(7) + C(8) + C(9) 
(^2))^/^ 

with a new denominator whose roots are much more manageable. In particular, the re-expansion leaves only the 
0{Ax^) terms, or those that are quadratic in the coordinate displacements, of the original denominator. From Eq. 



(12), we see that only s appears in the denominator, whose quadratic dependence on the coordinate displacement is 



simply 



r = g^^Ax^Ax^ + {usAx'^y + 0(Ax'^) (28) 



The second term is manifestly positive except at the location of the particle where it vanishes. The first term is 
not necessarily positive and may still potentially result in a vanishing denominator, in general. However, if in some 
coordinate system one chooses to associate the field point, x, with the particle position, x, so that they always share 
a common time coordinate (that is, t = t), then we have s^{t — t) — g^jAx^Ax^ + (ugAa;")^. Now, since gij is a 
purely spatial metric, its eigenvalues are all positive, so that the first term is unconditionally positive-definite and only 
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vanishes at the location of the particle. (See Appendix |B] for an explicit demonstration in the case of Schwarzschild 
coordinates). Thus, with a re-expansion of the denominator we achieve a simplification and, more importantly, we are 
also able to avoid the non-worldline singularities in the Haas-Poisson expression for the singular field given in Eq. ( 26 1 . 
The latter feature is essential for, say, a robust (3+1) application of the effective source approach. It is important to 
remember that to guarantee this, x and x must be on the same <-hypersurface, where t is the time coordinate in the 
specific coordinate system chosen to express $s- 



B. Periodicity of the singular field 

Although not strictly necessary, in spacetimes with axial symmetry it may be desirable to have an approximation 
to the singular field which is periodic in the azimuthal coordinate. There is no guarantee that that will be the case 



for the expansions (26) and (27); in fact there is not even any guarantee that <I>g(A0 = tt) = <i>g(A0 = — tt), i.e. that 
it is continuous across A0 ~ ±7r. 

Barack and Golbourn [JJ explicitly enforce periodicity by making the substitution Ai/)^ = 2(1 — cos A0) + 0{lS.(f)^). 
This was extended to higher order in Ref. j58j by making use of expressions involving cos(A0) and cos(2A0). However, 
both of these previous works only required replacements for even powers of A(/). In general odd powers of At/) can 
(and do) also appear. 

Among the infinitely many ways in which periodicity may be enforced for both odd and even powers, not all 
approaches are equal. For example, using replacements involving sinnAi/) proves to be a poor choice; since smn/S.(f) = 



at A(j) = ±7r for any integer n such replacements may lead to the denominator of (26) or (27) vanishing if Acj) = ±tt 



lies within the worldtube (Sec. IIICl) or the window function's region of support (Sec. IIIC2). Unfortunately, it is 
easy to see that no alternative choice for the functions used to replace odd powers of Acj) can avoid such extra zeros. 
That is, denoting the replacement for A0" by /„(A$), for any odd n the function /„ must have at least one zero 
somewhere in Aip € (0, 27r)|^ 

Given the two criteria: (i) $s(A0 — tt) — <bs{A<j) — —tt) and (ii) <I>s(A(/) — tt) is finite, we therefore propose 
a particular choice which satisfies both requirements and which has other practical advantages. We introduce the 
angular variables 

Q = sin i A0, R = sin Acj) (29) 

and rewrite even powers of Acj) in terms of Q and odd powers in terms of Q and R. This is easily achieved by 
expanding A(j) = 2arcsin(5 and A<j) — arcsini? for small Q and R, and making use of the identity R^ — 4(5^(1 — Q^) 
to give 

A0«i?(i + ^g-* + ^g5), A(/.2 « + 1q4^ A^b"" ^ mq^ (i + q^^ , a^^ « 16Q^ a^^^r^, (so) 

where terms of ©(A^^) and higher have been neglected. Not only does this replacement satisfy both criteria mentioned 
above, it also leads to relatively compact formulas for the partial derivatives 

= ^dg + (1 - 2Q^)dR, d^^ = ^(1 - Q^)dQQ - ^Qdg + ^(1 - '2Q')dQR - RdR + (1 - 2Q')^dRR (31) 

which appear in the wave operator (used when calculating the effective source). Moreover, this choice of variables has 



the subtle advantage of lending itself to minimal sensitivity to round-off effects close to the particle (see Sec. Ill C 3 
for an explanation of why this is important). For small At/), the substitutions of Refs. [43] and [55] are sensitive to 
numerical round-off, whereas this is not the case for our (Q, R) scheme. 



C. Calculation of the effective source 



With an approximation to the singular field at hand, we must now calculate a corresponding effective source. Before 
proceeding with the calculation, we will briefly mention some issues which one must be cognizant of. 



To see this, suppose that n is a (positive) odd integer. /„ must clearly satisfy the following properties (among others); 

1. /n is continuous 

2. fn{x) ~ x" for small \x\ (indeed, f{x) = x" + 0{x"^), where m is the highest power of x appearing in the expansion of the singular 
field) 

3. fnix + 2kiT) = fn{x) for any integer k 

Property [5] implies that if n is odd, /n > for small positive At/), and /„ < for small negative At/). Property |3] then implies that 
/„ < for At/) slightly less than 27r. Property [l] and the intermediate value theorem then imply that /„ must have a zero somewhere 
between At/) = and Atp = 2tt. 
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In general, the calculation of an effective source requires the computation of derivatives of $3 (2;)- When calculating 
these derivatives, one generally needs to be careful to take account of the fact that x' and x" vary with x since they 
must remain linked by a null geodesic [SS]. Fortunately, in the approximation $g of Eq. ( 12 ), by writing everything in 
terms of CTs, this dependence is made explicit in terms of the arbitrary point x which does not depend on x. However, 
since this dependence is only given as an expansion in e, it is an approximation which is only strictly valid in the 
limit e 0. For example, it is possible (and likely) that ( [l2| differs from the 'true' singular field at 0{e'\ yet in 
the limit e — > they agree. Similarly, the corresponding effective source has the correct value (i.e. 0) at e = 0, but 
contains differences from the 'true' effective source at 0(t). One must be careful to account for this when computing 
an effective source. 

In computing a covariant approximation to the effective source in Sec. |II C[ we made use of identities such as 
(T"(Ta = 2cr and <J°^(Ja = 2(7. Furthermore, we re-expressed higher derivatives of u in terms of their covariant expansion 
about X. However, once coordinate expansions are introduced, these identities and expansions are no longer exact - 
they are only valid up to the order of the coordinate expansion. In computing the singular field, this is not an issue 
since we are only interested in computing the value of the self-force at the particle, in which case the errors vanish. 
Unfortunately, the effective source is required not just at the particle, but also in a region surrounding the particle 
where the errors are no longer zero. It is therefore not possible to make use of these simplifications when calculating 
a coordinate effective source (at least not without taking care that $g and its derivative evaluated at the particle are 
unchanged). 

With these issues in mind, there are now two choices on how to proceed with computing S'off = — D^g- We will 
investigate both of these in turn in the following sections. 



1. World-tube method 



Barack and Golbourn [44j propose a precise method for computing an effective source. They introduce a world-tube 
around the particle. Inside the world-tube one solves for $r and outside one solves for <I>rct, which is now a solution of 
the homogeneous wave equation in this region. By imposing the boundary condition $rct = "I'r + one can ensure 
that the system as a whole is consistent. 

They look for a 'puncture' field - an approximation to the singular field which depends only on the spatial position 
of the field point, with all time dependence encapsulated in the particle motion, 

<i>p{x\x^t),u''{t)). (32) 

In doing so, they effectively t — t, i.e. fixing x to depend on x in the sense that their time coordinates are equal. 
Recall that x is arbitrary and does not have any required dependence on x. Their choice is therefore valid and 
consistent with the singular field computed in Sec. |III A| In particular, their choice of puncture function 

*p = (33) 
^y {gij + UiUj)Ax'-Ax^ 

corresponds exactly to the first order singular field given here (with t — t). In Ref. [BS], Barack, Golbourn and Sago 
proposed an improved puncture function, which again is equivalent to the second order singular field given here. 



Computing higher order puncture functions is straightforward: one takes the expansions (26) or (27) at the desired 



order and sets At = 0. For example, a fourth order puncture function for a particle in circular equatorial geodesic 



motion around a Kerr black hole is given explicitly in Sec. IV B Dolan and Barack ,47] have recently made use of a 



similar fourth-order puncture computed in this way to calculate the self-force on a particle in a circular geodesic orbit 
about a Schwarzschild black hole and Dolan, Barack and Wardell [S5] extended this to the Kerr case. 

Given this puncture field, the computation of an associated effective source is straightforward. One simply calculates 
an expression for the wave operator in the coordinates in which $p is given and applies this wave operator to , 
noting that spatial derivatives act only on x', while time derivatives act only on x^{t) and u°'{t). 



2. Window function method 



In a numerical 3-1-1 evolution code, it is most straightforward to solve for $r everywhere, requiring SoS everywhere 
on a 3D spatial slice. This would be problematic wherever <I>s is either not defined or where its series expansion 
diverges. Vega and Detweiler |45) propose a solution which involves the introduction of a window function, W{r), 
which smoothly transitions from a value of 1 at the source to far away. In effect, one is then solving for near the 
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FIG. 3. Coordinate series expansion approximation to the effective source close to the particle. The exact effective source 
(solid blue line) is well approximated in the region Ar < 0.05M by the first term in its series expansion (dashed purple line). 
Including the second (dashed gold line) and third (dashed green line) terms further improves the agreement for points farther 
from the particle, with the curves being almost indistinguishable from the exact curve. 



particle and for <i>ret far from the particle. For this to provide a consistent C'^ source, we must impose a restriction 
on W: at the particle it must be 1 and at least its first three derivatives must be 0. Having introduced this window 
function, the effective source is then given by 

S-eff = V{W^s)- (34) 

In |46| . Vega et al. use an expression for ^efi which is (i.e. continuous but not differentiable) , limiting the convergence 
of their finite differencing scheme despite their use of 8-th order spatial differencing. A smoother source would be 
advantageous in that it would give a higher convergence order without the need to construct a more complicated 
finite differencing scheme to deal with the non-smoothness of the source. In this work was extended to include 
the back reaction from the self-force into the evolution. This latter calculation made use of the re-expanded singular 
field described above, evaluated in Kerr-Schild coordinates with the choice ^ks = ^ks- 



3. Evaluation of the effective source very close to the particle 



Severe round-off errors may be incurred when evaluating the effective source very close to the particle. Applying 
the wave operator to the singular field results in many terms that scale as 0{e~^), which evaluate to large quantities 
as e — 0. However, we show in the analysis of Sec. |II C| that, at the order of our present approximation to the singular 
field, all of these terms cancel to leave an over-all effective source that scales as 0(e). As was already pointed out in 
[42], this is a prototypical example of catastrophic cancellation that is often encountered in numerical work. There 
are two solutions to this problem which have been found to work well. We describe each approach in detail below and 
note that the choice of which scheme to use is dependent on the problem at hand. For simpler configurations (with 
more manageable expressions) a series approximation may be appropriate, whereas in cases where more unwieldy 
expressions appear it may be more straightforward to use numerical interpolation. 

Series approximation close to the particle 

Given that cancellation is only an issue for points very close to the particle (typically at a distance of < 0.05M), it 
is reasonable to replace the full effective source in this region with an approximation which is valid for points a small 
distance from the particle. In particular, by replacing the full effective source with its series expansiorj^ one obtains 
an expression for the effective source of the form 

A _ C^(12) + (^(13) + ^(14) , 4 /o.^ 

5eff ^^^^^72 + '^(^ )■ (35) 

This expression is manifestly C(e), with all divergent terms having been cancelled analytically. For the small region 
where catastrophic cancellation arises, it is sufficient to take only the first term, j^i^^njs, which is 0{t). The inclusion 
of subsequent terms would only be necessary if an approximation was needed in a much larger region (see Fig. |3| . 
There is a potential disadvantage to this scheme, however, in that it involves the evaluation of the twelfth order 



Although our approximation to the singular field is already written as a truncated series expansion, the effective source which is computed 
by applying the wave operator to it is not. For example, as can be seen from Eqs. | |45| l and ( |53| there are several terms which depend 
on the location where the source is being evaluated and which are not written explicitly as series expansions. 
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polynomial, d(i2), which may be quite computationally expensive. Fortunately, since this is to be implemented only 
in a very small region around the particle, the overall computational burden this adds is likely to be minimal. 
Interpolation close to the particle 

Another solution to the catastrophic cancellation problem is to compute the effective source via interpolation. We 
take advantage of two important facts: (1) SeS = on the worldline 7, and (2) SeS is smooth everywhere except on 
7, where it is just C° . The first fact gives us an exact data point for the effective source, while the second justifies 
an assumption that interpolation might be sufficient. We identify a small region TZ around the particle location 
outside of which the effective source is computed reasonably well. If the effective source is required inside TZ, say at 
then it is first evaluated at selected points along a "coordinate ray" outside this region. Using these values and 
Scsix"^ = a;*) = 0, where a;* are the spatial cooordinates of the location of the particle, we then interpolate to y*. (All 
coordinates here are purely spatial in compliance with the restriction mentioned in Sec. [Ill A[ when evaluating the 
effective source, all field points must be on the same t-hypersurface as the particle.) 

More concretely, consider 5(A) := S'off(a;*(A)) as a function of A, along the coordinate ray given by x'(A) = 

+ A(y' — x'). If e TZ, then to compute Scsiy^) = S{X — 1), we interpolate using a few evaluations of S{Xj) [such 
that x^lXj) i TZ] and ^(A = 0) = 0. 

Obviously, there is considerable freedom in how to implement a specific interpolation scheme and in what to choose 
for the size of the interpolation region TZ. The results reported in [43] appear to be very robust with respect to the 
various choices we have tried. 



D. Specific scliemes 



There are three commonly applied approaches which may be used for solving the wave equation for the regularized 
field. These methods solve for the regularized field in 1 + 1, 2 + 1 and 3+1 dimensions, eliminating the other dimensions 
through a decomposition in suitable basis functions. In the 1 + ID approach, the spherical harmonic basis is chosen 
and the decomposition is done into I and m modes, while in 2 + 1 dimensions the decomposition is only done into m 
modes. There is a trade off between having to evolve a field in higher dimension (and all the difficulties of poor scaling 
with resolution that goes with it) and requiring the calculation of a large number of modes. It remains to be seen 
which is the better choice; the conclusion will most likely depend on the particular configuration under consideration. 
Nevertheless, the calculation of the singular field and effective source proceeds in the same way. Both are calculated 
as 3 + 1 dimensional quantities using the methods described in the previous sections. In the 2 + ID m-mode scheme, 
they are then decomposed into m modes by performing an integration over the azimuthal coordinate, <p: 

^s{^r,A0,Acj),At)e-'^Ucl) ^cff (m) - / S,s{Ar, A9, Acf>, At)e-'^^dcp. (36) 

-TT J —TT 

For the 1 + ID m-mode scheme, a second integration is performed over the inchnation angle 9: 

4's(i,m)=/ / ^s{^r,A9,A,p,At)Y*^{9,^)d9d(^ 

J-TT Jo 

/TT nTT 
/ S,s{Ar,A9,A(l>,At)Yi:^{9,cl))d9dcb. (37) 

In practice it may be most straightforward to do the integration numerically. As a result, the calculation of the 
singular field and effective source may dominate the runtime of a 1 + ID or 2 + ID code. 



IV. EXAMPLES 



In this section, we give examples of the singular field and effective source in some specific scenarios. We consider 
in detail the case of a scalar charge undergoing circular, equatorial, geodesic motion in both Schwarzschild and Kerr 
spacetimes. Note, however, that the methods developed here do not depend on the symmetries present in these 
configurations. They are equally effective in other spacetimes and for generic geodesic motion. For these more generic 
configurations, the results are most easily given in electronic form. For this reason, we provide expressions for more 
generic configurations online |60j and give here only explicit examples for simple configurations along with plots for 
more generic configurations. 
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A. Circular geodesic in Schwarzschild 

Given the Schwarzschild metric in standard coordinates, 



ds^ = - 1^1 j dr +1^1 j dr^ + r^dr + sin^ (38) 

we foUow the prescription of Sec. |III| to obtain a fourth order approximation to the singular field of the kind given 
in (27). In general, this will be a function of the field point, x = {r,9,(j3,t), the world-line point x — {f,9,(j),t) and 
the particle four- velocity m° = {u^ ,u'^ jU*^). We may use the spherical symmetry of the spacetime to enforce that 
the motion lies in the equatorial plane, i.e. 6 = 7r/2, = 0. In order to obtain sufficiently compact expressions to be 
given here, we make the further assumption that the motion is circular, i.e. [70j 



f = constant w'' ==0 it'^ = -\ — — — u* = W . (39) 

Finally, we use the freedom in the choice of i to set the field point and world-line point to be at the same coordinate 
time, i.e. 

t = t 4) = nt, (40) 

where 

is the orbital frequency. Combining everything, we obtain a fourth order approximation to the singular field of a 
scalar charge on a circular, equatorial orbit around a Schwarszchild black hole: 



where Ar — r — f, A6 — 6 — 7r/2, Q — sin (^(0 — f2i)) and where the non-zero coefficients, Uijk and bijk are functions 
of the orbital radius, f, and are given by 

64f6(2M-f)3 _m^{f~2M)^ _12f^{f~-2M) _ _g _48f^{f-2M) 

(r--3M)3 ' ""^4- (r--3M)2 ' ""'^^ " f - 3M ' "o^o - r , «204 - (-.3^)2 - 

24f5 _ 3f5 _ 12f4 _ 3f4 _ f3 

^^^^-"3M^' ''^^°-"2M^' 6M2-5Mf + r-2' '^^^^ " (r - 2M)2 ' '^s"" " (r - 2A/)3 ' 

32f5(M - f)(f - 2M)2 _ 8f5 (8M2 - lOMf + 3f^) _ 2f5(3f - 5Af) _ 

(3A/-f)3 ' "12^- (f-3Af)2 ' "1^^- 3Af-r- ' "i^o - -y, 

_8f4(3A//-2f) _ 8f4 _f4(5Af-2f) _ 2f3 _ f^^{f-4:M) 

(f-3A/)2 ' "322-3Arr^' ^340- 2(f^_2A/)2 ' «502 - - _ ^^^^^ , «520 - -^^^r^^, 

Arf2 32Afr5(2M - f)3 16f5(f - 2Ar)2 (SM^ - 4Mf -f r^^ 

^700 = m,T\d ' '^008 — TTTTT 3:^ 7 ^^026 ^ 



2(f-2Ar)4' '^^^^ (3Ar-f)3 ' (3A//-f)3 

2f5 (222Af4 - 459Af 3f 346M^f^ - 112Aff3 -)- 13f'') _ (-30Af3 + STAfV - 29A//f2 + 4f3) 

3(f-3A/)3 ' " 3(f - 3 A/) 2 



(6M2 - 9Mf r^) 8f4 (-30Af 3 + 35A/2f - 16Mf^ + 3f 



3\ 



a080 — -r. „ J ?r73 ; ^206 



72M-24f ' (f-3Af)3 
2f4 (37A/2 - 40Arf + 13f2) _ (294A/3 - 498Af + 259Mf^ - 41r 



0-224 — 7^ „ , 1 ^242 



(f-3Ar)2 ' 6(2A/-f)(f-3Af)2 

(48Af2 - 57Aff + llf2) 2f3 (-65Af3 -|- 74M2r - 26Arf2 + 3f3) 



'^260 — nA /g 7i/r9 r 71 J- I -9\ ' ''404 



24(6Ar2-5A//f + f2) ' ™* (2A/-f)(3M-f)3 
_ -50M3f3 _^ 72M2f4 _ 30Mf5 4f6 _ 5f3 (-9^2 _ ^oMf -f 2r2) 

(6M2-5Mr + r-2)2 ' " 24(r- - 2Af)3 
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Mf2 (31M2 - 37Mf + 8f2) _ Mf^ (19^2 - 21Mf + 4f2) 



^^602 — o,,rN9 ' "620 



2(2M-f)3(f-3M)2 ' "-^^ 8(3M-f)(f-2M)4 
M2f(f - M) 16Mf4(f - 2M)2 (29M3 - 25M2f + 3Aff 



"800 — ./o7i/r -\.;/o7i J =V' "108 



4(2M-f)5(3M-f)' (3M-f)5 

4f4(2M - f) (97M4 - 86M3f + 27M2f2 - SMf^ ^_ 2f^4-) 
(f- 3M)4 



(312M4 - 351M3f + 193M2f2 - 73Mf^ + ISf") {5AM^ - bUPf + 31Mf^ - 8 



ai44 — 777^ „ , ,r\-i 7 "162 — 



3(f-3M)3 ' 12(r-3M) 

f5(3M + f) 4f3 (139^4 - l63A/3f + 78Af2f2 _ 27Mf^ + 5f 



2 



4^ 



"I8O — ,n/o , ,r "306 



48(3M - r) ' (f - 3Af )4 

f3 (-357Af 4 + 434Af3f - 230Af2f2 + 78A4"f3 - 13f4) 
0324 — 



(2Af-f)(3A/-f)3 

f3 (732A//4 - 1074Af3f + 677Af2f2 - 239Af + 38^"*) _ {-66M^ + 69M^f - iSMf^ + Uf^) 

" 12(6Af2_5Aff + f2)2 ' " 48(3Af-f)(r-2Af)2 

_ Aff2 (l95Af3 - 207A/2f + 83Aff2 - 9f^) _ (-279Af4 + 307Af¥ - 86Af2f2 _ i2Aff3 _^ 4^4-) 

" (3A/-f)3(f-2A//)2 ' "^^^ " 4(f-3A//)2(f-2A//)3 ' 

2 (132A4^4 _ 75M3f - 41Af2f2 _^ 32Mr3 „ 2f4^ (89A4^3 _ 73^2^ + lOAf f2 + 4f3) 



r 

"540 — ^ -o^oTi J =w"^ r.7. ^' "702 



and 



48(3Af-f)(f-2Af)4 ' '"^ 4(f-3Af)2(f-2Af)4 

Mr (- 19Af 3 + 7Af2f + 6Af f 2 - 4f 3) Af 2 (2Af 2 _ 2Aff + f^) 

" 16(2Af-f)5(3Af-f) ' " 8(3Af-f)(f-2Af)6 ^^^^ 

4f2(f-2Af) , -2 , ^ 

^002 = — z — iTji—, O020 = r , &200 = z — ^ (44) 
r — 3Ai r — 2Af 

Next, we compute the effective source corresponding to this singular field. The wave operator in Schwarzschild 
coordinates is given by 

_ / r \ 92 fr-2M\ 2{r - M) d 1 Id Id 

^'^^^ ^^\r-2M)W^\ r J ^ ^ d^^V^W^ r2 tan(6') M ^ sin^{e)W ^ ^' 



- , . . . (46) 



Applying this to (42), we obtain an effective source of the form 

.(4) _ /(Ar,Ag,g) 

'eff — , 

(E::^f.t'o-'&..^Ar^A0^AQ'= 

where /(Ar, A0, Q) is a polynomial in Ar, A0, Q (and contains terms involving tanA^? and secA0) and the are 
the same as those in the singular field. 

In Fig. [4] we plot the first, second, third and fourth order singular field and corresponding effective source for the 
case of a particle in a circular orbit at f — lOM in Schwarzschild. All four cases have a visually similar singular field. 
This is not surprising given they share the same singular behaviour and only differ in higher order corrections. The 
corresponding effective source, however, is very different. As expected from the discussion of Sec. |II C[ at first and 
second order the effective source diverges at the particle, while at third and fourth order it is finite. Figure [5] shows a 
zoomed in view of the effective source in each case, along with a slice along the radial direction, passing through the 
particle. From this we see more clearly the behaviour of the effective source near the particle: at first order it is C~^, 
at second order it is C~2, at third order it is C^^ and at fourth order it is C". 

B. Circular geodesic in Kerr spacetime 

To compute the singular field and effective source in Kerr spacetime, we consider its metric in Boyer-Lindquist 
coordinates, 

/ 2Afr\ ,0 4aAfrsin2 6l , ,, S,, ^,,0 2Afr(r2 + a2)\ . „ , „ 
ds^ = - [1 — j df dtdcj) + —dr^ + Y^dO^ + ( A + ^ ) sin^ ed(j)^ (47) 
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FIG. 4. Singular field (top) and effective source (bottom) along the equatorial plane for a particle in a circular orbit around a 
Schwarzschild black hole. From left to right: first, second, third and fourth order cases are shown. Note that in the third and 
fourth order cases, we used the method described in Sec. [Ill B] to ensure periodicity in 0. 




FIG. 5. Close-up view of the effective source along the equatorial plane (top) and along a radial slice through the particle 
(bottom) for a particle in a r = lOM circular orbit around a Schwarzschild black hole. From left to right; first, second, third 
and fourth order cases are shown. 



where 



S = + cos^ 9, A = - 2Mr + 



(48) 



As in the Schwarzschild case, in order to obtain sufficiently compact expressions to be given here, we assume that the 
motion follows a circular, prograde equatorial geodesic, i.e. |70j 



u"- = 0, = 0, u"^ 



r = constant. 



TT 

2 



aAI + VMr^ 



V?-^ - 3Afr + 2a^/Mr 



(49) 



MrJir"^ -3Mr + 2aVMr) 



We also use the freedom in the choice of t to set the field point and world-line point to be at the same coordinate 
time, i.e. 



t = t 



= nt 



where 



M 



aM + VMf^ 



(50) 
(51) 
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FIG. 6. Particle following a circular equatorial geodesic around a Kerr black hole with spin a — 0.99M. Left to right: (1) 
fourth order singular field, (2) fourth order effective source, (3) close-up view of the effective source and (4) effective source 
along a radial slice through the particle. In (4), we also show the Schwarzschild result as a dashed red line for comparison. 
Note that we used the method described in Sec. |III Bj to ensure periodicity in 0. 



is the orbital frequency. Combining everything, we obtain a fourth order approximation to the singular field of a 
scalar charge on a circular equatorial orbit around a Kerr black hole: 



7/2 



(52) 



where Ar = r — f , = 9 — 7r/2, Q = sin (|(0 — ^t)) and where the non-zero coefficients, a^-fc and fci^fe are functions 
of the orbital radius, f, and the spin parameter, a, and are given by taking the expressions in Ref. ^58j, making the 
change of variables A0 Q and re-expanding as described in Sec. |III A| 

Next, we compute the effective source corresponding to this singular field. The wave operator in Kerr (Boyer- 
Lindquist) coordinates is given by 



□ 



BL 



2Mr(a'^ 



(a2 + r2 - 2Mr) (r"^ + a? cos^ 9] 
1 cot6' d 



at2 



2Mr) 92 



(r2 



r2 + a? cos2 9 dr^ 
2Mr + a2 cos2 9) csc^ 



2{r ~M) d 
2 + a2 cos2 9 dr 

92 



r2 + a? cos2 9 d9^ + cos^ 9 89 {a^ + r2 - 2Mr) (r2 + o? cos^ 9) dcj? 



AaMr 
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(a2 -t- r2 - 2Mr) {r^ + cos^ 9) d<j)dt ' 



Applying this to (52), we obtain an effective source of the form 

/(Ar,A0,g) 



q(4) 
'cff 



11/2 ■ 



(53) 



(54) 



where /(Ar, A6', Q) is a polynomial in Ar, A6', Q (and contains terms involving tanA0 and secA0) and the bijk are 
the same as those in the singular field. 

In Fig. |6] we plot the fourth order singular field and corresponding effective source for the case of a particle 
in a circu lar orbit at f = lOAf around a Kerr black hole with spin a = 0.99Af . As expected from the discussion of 
Sec. II C[ the fourth order effective source is finite and continuous, i.e. C". In the rightmost figure, we compare against 
the equivalent case in Schwarzschild. Both cases are qualitatively remarkably similar, only differing significantly in 
magnitude close to the black hole. 



C. Generic geodesic in Kerr spacetime 

To illustrate the power of the method developed here, we now consider a more generic configuration. We choose an 
arbitrary timelike geodesic of the Kerr spacetime and compute the singular field and effective source at a point along 
that geodesic. In particular, we make the choice 

f=10M, a = 0.99M, M = l 

2 

VMt g 1 ^ 



r 



\/r2 - 2>Mr + 2a\/Mr' 2 



= u'^ = Mm* (55) 
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FIG. 7. Particle following generic geodesic around a Kerr black hole with spin a = 0.99M. Left to right: (1) particle's world-line 
(solid line) with position at which the singular field and effective source are computed indicated by a black dot, (2) fourth 
order singular field, (3) fourth order effective source, (4) close-up view of the effective source. Note that we used the method 
described in Sec. [Ill fj] to ensure periodicity in 0. 



with u* being determined by the normahzation of the four- velocity, u^u" = —1. The computation of the singular 
field and effective source proceeds exactly as in the circular orbit case, the only difference being that the resulting 
expressions are larger. In fact, they are too large to be useful in printed form. Since they are relatively manageable 
with computer algebra, however, we have made them available as Mathematica code |60j . 

In Fig. [7| we illustrate the behaviour of the singular field and effective source for this configuration. The leftmost 
plot shows the geodesic over several orbits, indicating that it is both inclined and eccentric. The black dot on this 
plot indicates the point f = lOAf, 6 = 7r/2, ^ = at which the singular field and effective source in the subsequent 
plots is computed. The second plot shows the singular field along the equatorial plane. The third and fourth plots 
show the effective source along the equatorial plane. As expected, this fourth order effective source is continuous, but 
not differentiable at the particle. 



V. DISCUSSION AND SUMMARY 



In this paper, we have developed an approximation to the singular field of a point scalar charge to quadratic order 
in the distance from the charge. This is sufficient to give second order convergence in the grid spacing for 3 -I- ID 
numerical calculations and to give m~'^ and convergence in the m-mode and ^,m-mode schemes, respectively. To 
go to higher order (for better convergence) one would need to: 

1. Calculate the higher order terms in the coordinate expansions of da- This is a recursive calculation and the 
expressions get more unwieldy as the order increases. However, the calculation method is general and only 
limited by computational power. 

2. Calculate higher order corrections to ctq'u" and (Ja"U°' . This is straightforward using higher order covariant 
expansions of a and its derivatives [57j, which are easily obtained to much higher order than is needed here 
using non-recursive methods |71j . 

3. Calculate higher order terms in the series expansions of t^(x, x'), /7(x, x") and T^(x, ^(t)). These are also easily 
obtained from the semi-recursive methods of Ref. [71] . 

The calculation of a higher order singular field and effective source is therefore a straightforward (if somewhat tedious) 
process. With the expressions becoming more unwieldy at each order, one must balance the calculation effort against 
the benefits of doing so. It seems likely that the fourth order approximation presented here is the 'sweet spot', giving 
reasonably good convergence with modest computational difficulty. 



As shown in Sec. (IV C), this method works for very general motion in the Kerr spacetime. Furthermore, although 



no explicit calculations have been done here for other spacetimes, it is clear that Eq. ( 12 1 is valid in any spacetime. It 
would therefore be straightforward to apply this method to any (not necessarily Ricci-fiat) spacetime. As the primary 
motivation of this work has been to study rotating black holes, however, we have chosen to only consider Kerr and 
Schwarzschild spacetimes in detail in this work. 

The methods presented here are useful for computing expressions for the singular field and effective source for 
generic configurations. The actual evolution of a wave equation with this source, along with the calculation of the 
self-force should be explored separately. Here, we simply note that we have implemented two separate numerical 
evolution codes using the singular field and effective source presented here: one uses the window-function approach 
with a 3 -|- ID numerical evolution; the other uses the world-tube approach followed by an m-mode decomposition and 
a separate 2-t-lD numerical evolution for each m. We have verified that both codes give correct results (as determined 
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by comparison with frequency-domain calculations). Further details of these codes will be presented elsewhere, with 
some results already having been published. Using a separate numerical code, Dolan an d Bara ck [47] evolved the 



2 + ID scalar wave equation (with the singular field and effective source as given in Sec. TV A ) for a particle in a 
circular orbit around a Schwarzschild black hole. This calculation was subsequently extended to the case of circular 
orbits in Kerr spacetime in [581 172| with further progress toward generic configurations in Kerr spacetime under way. 
In a recent work [43| , the effective source presented here was used to self-consistently evolve the orbit of a point scalar 
charge in the Schwarzschild spacetime, incorporating the back reaction from the self-force into the evolution. 

One major issue remaining in the effective source approach is the computational efficiency of the source calculation. 
For the approach to be of practical use, its calculation must be sufficiently fast that it does not have a prohibitive 
impact on the run time of a numerical code. This is a serious concern - the expression for the fourth order effective 
source may be dramatically larger than a finite difference representation of the wave equation, for example. Some 
steps have been taken in this paper to improve the efficiency of the source calculation. In Appendix [Cl we discuss 
some specific methods for evaluating the effective source as efficiently as possible. As mentioned in Sec. |III A[ we have 
also made use of specific choices for the singular field in an effort to minimize the size of the resulting expressions. 
Despite these efforts, the reality is that the calculation of the effective source will considerably affect the run time of 
a numerical code. 

Fortunately, there remain several possibilities for further optimization. The advent of GPU (Graphics Processing 
Unit) computing has allowed for dramatic performance improvements in certain applications. It seems likely that 
the embarrassingly parallel nature of the effective source calculation on a grid of points is an ideal candidate for 
implementation in a GPU programming framework such as CUDA or OpenCL. Given other applications have seen 
speed-ups by 1 to 2 orders of magnitude [73] , it is not unreasonable to expect similar performance gains for effective 
source calculations. 

There is yet another intriguing prospect for improving calculations involving an effective source. As discussed in 
Sec. |II C 4 the effective source may be viewed as merely a correction for the fact that the singular field is not known 



exactly. This begs the question of whether the singular field could be calculated exactly on a world-tube boundary. 
Not only would this improve convergence in a numerical code (arbitrarily high convergence in grid spacing, exponential 
convergence in I or m mode sums), but it would also negate the need to calculate an effective source at all. The 
entire computational cost of implementing the effective source approach would be in the computation of the value of 
the singular field on the boundary. While an exact calculation of the singular field may not be realistic, one should 
recall that from a numerical perspective a value which is correct in the first 16 digits of a double precision number is 
effectively 'exact' in that further refinements do not change the result. Given the availability of high order expansions 
of the Green function [S31 [7TJ [7iti75] along with the fact that multi-domain spectral methods [29, 30] or adaptive 
mesh refinement [41] [79] allow the world-tube boundary to be placed very close to the particle, it seems like this may 
be a plausible approach, although further investigation is required to determine whether this is truly the case. 



Yet another potential optimisation arises from the covariant treatment of Sec. II C Near the particle, the covariantly 
re-expanded effective source is a reasonable approximation to the 'correct' effective source. However, given that it only 
requires a first order coordinate expansion, it is dramatically more efficient to evaluate numerically. Furthermore, 
as the divergences are cancelled analytically, it effectively avoids any need for concern about delicate numerical 
cancellations. Lastly, as the singular field is constructed in such a way that it and its derivative (i.e. the self-force) 
evaluated at the particle are insensitive to these covariant re-expansions, it is plausible that using the covariant re- 
expansion throughout the world-tube may be possible. This may lead to an 'incorrect' regularized field away from 
the particle, but with sufficient care could potentially still give the 'correct' value for the field and its derivative at 
the particle. 

This work focused on the case of a scalar charge moving in a background spacetime. Of arguably much more 
interest are the cases of gravitational or electromagnetic charges. Fortunately, the calculation strategy remains largely 
unchanged. One can make use of an analogous Detweiler- Whiting gravitational or electromagnetic Green function 
which has the same Hadamard-type structure. It will still include functions U{x, x')^ b' and V{x, x')^b' which are 
analogous to their scalar variants and may be calculated in the exact same way [71| . Furthermore, the world function, 
a, and its derivatives will remain unchanged from the scalar case. The full details of this calculation will be developed 
in a future work. 



In fact, the singular field used by Dolan and Barack differs slightly from that of Sec. |IV A| Nonetheless, it was computed using the same 
methods and differs only at higher order than the order of the approximation. 
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Appendix A: Covariant expansions 

In this appendix, we develop covariant expansion expressions for the biscalars U{x,x'), U{x,x"), (Ta'{x,x'), 
<7a"{x, x') and /J V{x, z{T))dT appearing in Eg. ([TT| ). We eventually seek expansions about the point x. In doing so, 
we follow the strategy of Haas and Poisson [371 155] : 

• For the generic biscalar A{x, z{t)), write it as A{t) = A{x, z{t)). 

• Compute the expansion about t ~ f. This takes the form 

A{t) - A{f) + i(f )(t - f) + \A{f){T - f)2 + • • • , (Al) 

where A{f) = A-g^u^ , A{f) = A.^^^u^u^ , 

• Compute the covariant expansions of the coefficients A(f), A{f)^ ■ ■ ■ about f. 

• Evaluate the expansion at the desired point, e.g. A{x') = A{x,x'). 

• The resulting expansion depends on r through the powers of r — f . Replace these by their expansion in e (about 
x), the distance between x and the world- line. 

A key ingredient of this calculation is the expansion of A = r — f in e. This expansion was developed by Haas and 
Poisson [37] to sufficient order for the present calculation for the particular choices = v — f and A_ = u — f. 
They found 

A± = (r ± S) T ^'^-^ Ruaua T !^ [(r ± s)R^aua\u " Ruaua\a\ + 0{t'). (A2) 



1. Expansion of U{x,x') and U{x,x") 



We now compute expansions of U{x,x') and U{x,x") about x. Both calculations proceed in the same way and 
require the expansion of U{x, x) about x which is given by |52J: 

U{x, x) = ^^l\x, x)^l + ^R,„ ~ ^Raa\a + O(e'). (A3) 
Writing U{t) = U{x, z{t)), where t stands for either u or v, we compute its expansion about t ~ f: 

U{t) = U{f) + t7(f)(r - f) + \u{f){T - ff + \u{f){T - ff + 0{e^), (A4) 

2 D 



where 



U{f) = 1 + ^i?.. - + O(e^) (A5) 

U{f ) = U,o.u'' = ^Rua - ^Rua\a + ^^..|u + O(e') (A6) 

U{f) = U.^pU"u^ = ii?„„ - ^Ruu\a + lRua\u + 0{e^) (A7) 

jjif) = a^-p^u^u^u^ = + 0(e). (A8) 
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Substituting Eqs. (A2) and (A5|-(A8| into (A4) and evaluating at t = {u,v}, we get our final expression for the 
expansion of C/_ = U{x, x') and U+ = U (x, x") about x: 



U±=l 



1 

12 



1 

24 



-R. 



l„„ + 2(7 ± s)Ru„ + (r ± sfRm 

va\cr + {Racr\u ~ '^Rua\a){^ ± s) + {2R^„\^ — i?„„|o-)(r ± s)^ + -R„„|„(r ± s)^ 



(A9) 



I \ ^ I • 

The first term here is 0(1), the second term is O(e^) and the third term is 0{e^). Note that for vacuum spacetimes, 
these become U{x,x') = 1 + 0{t^) = U{x,x"), as is to be expected. Additionally, note that the difference between 
U{x,x') and U{x,x") first becomes apparent at C(e^). 



2. Expansion of a^'u" and (Ta"U° 

Haas and Poisson give expansions for aa'u" and aa"U°' . They are: 

r — s 



(T„'U — S - 



(Tr,"U — S 



-Ri, 



s 24 s 



-Ru 



6- - ^U(JU(J c\ A — 

s 24 s 



[(? - s) (r + 2 s) i?„<,„^|„ - (r + s) Ruaua\a\ + C(e^) 
[(? + s) (r - 2 s) - (r - s) Ruaua\a\ + C(e^), 



(AlO) 



(All) 



In these expressions, the first term is ©(e), the second is Oi^e") and the third is Oit^). Note that difference between 
crc'U" and (Ja"u'^ only becomes apparent at 0{e^). 

3. Expansion of j^V(x,z{T))dT 



The expansion of the tail term in Eq. (11) poses an additional potential difhculty because of the integration over 



a portion of the world-line. However, expanding V{x, z{t)) about i, the integration becomes a trivial integration of 
powers of t. 

In the following, we make use of the expansion of V{x, x) about x, 



v{x,x) = i(e- l)i?- ^(e- i)i?|, + 0(6^). 



(A12) 



We now proceed, as before, by defining V{t) = V{x, z{t), where r lies between u and v, and computing the expansion 
about T — f: 



where 



V{t) = V{f) + V{f){T-f), 



V{f) = l{^-l)R-\{^-l)Rl, + 0{e') 
y(f) = ^„u" = i(e-^)i?|„+0(e). 



(A13) 

(A14) 
(A15) 



The integration along the world-line is now straightforward since the only dependence of the integrand on r comes 
through the factor r — f. Performing the integration and substituting Eqs. (A2), (A14) and (A15) into the result, we 
get our final expression for the expansion of V{x, z(T))dr about x: 



V{x, z{T))dr ^{i-\)Rs+\{i- \){R\u-rs- R\a s) + Oie"). 



(A16) 



The first term here is 0{e) and the second term is ©(e^). Note that for vacuum spacetimes, V{x,z{t)) = 0{e^) and 
this term does not contribute to the singular field until 0{e^). 
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Appendix B: Leading-order piece of the coordinate expression for 

As we discuss in Sec. |III A[ it is highly desirable for a coordinate representation of the effective source to not diverge 
anywhere. Unfortunately, a common feature of series expansions is that they have a finite region of validity (for Taylor 
series, this is denoted by their radius of convergence). Outside this region, spurious singularities tend to appear. As 
already indicated, a re-expansion of the denominator of the singular field (leaving only its quadratic leading-order 
dependence on the coordinate separation and bringing all higher-order terms up to the numerator) allows one to avoid 
most potential singularities away from the position of the particle. Here we demonstrate this explicitly for the case 
of Schwarzschild coordinates. 

In these coordinates, the leading-order dependence of the denominator (essentially given by = (g"^ + u"u^)(TaO'j;s) 
on the coordinate separations is 

s2 = - Fil - F{u*f){Atf - 2u*M'^(Ai)(Ar) - 2FR'^u'u'f'{At){A(t)) 

+ R'iAef + ^^^^^^{Arf + '^u^u^{Ar){Acj>) + R\1 + R^u'^f^iA^f, (Bl) 

where At — t — t, Ar = r — R, etc. (recalling that barred coordinates refer to the position of the particle), R is the 
radial position of the particle in Schwarzschild coordinates, and _F := 1 — 2M/R. 
If we take At — 0, this reduces to 



s2(A< = 0) =i?2(Ae»)2 



(Ar)2 + -^u'-u'^(Ar)(A(/.) + R^{1 + R^{u'^f){A(t>) 



p2 y ' • p 



(B2) 



The condition At — imposes that the position and four- velocity of the particle are evaluated at the same coordinate 
time as where the effective source is evaluated, or in other words, the particle location and field point need to be at 
the same t-hypersurface. 

All except the cross term cx (Ar)(A0) are manifestly positive-definite. The combination of terms in the square 
brackets, however, can also be shown to be positive-definite; it is a quadratic form in {Ar, A(p}: 

A(Ar)2 -I- B{Ar) ( A0) + C( A^)^ (B3) 

where 

^^(J^+KT)^ ^^2^^,^^^ C = i?2(l + i?2(^0)2) (B4) 

This will be positive-definite if — 4AC < 0. The condition is easily verified to reduce to 

0<F+{u'y + FR'^{u'l'f, (B5) 

which is true for F > 0. Thus, s^(At = 0) is positive everywhere except that it vanishes at the position of the 
particle. This implies that the re-expanded singular field, which keeps only the quadratic dependence on the coordinate 
separation in its denominator, diverges only at the location of the particle, and consequently, that the corresponding 
effective source is regular everywhere else. 



Appendix C: Efficient numerical computation of the singular field and effective source 

The calculation of numerical values for $s ^^nd S'cff requires the numerical evaluation of their coordinate expansion. 
This amounts to numerically evaluating a multivariate polynomial (in {Ar, A6, A(f>)) with coefficients which are 
potentially complicated functions of the particle's location and four-velocity. Furthermore, in a numerical code this 
must be done at every point on a 3D gridj^ Clearly, it is crucial to make this evaluation as efficient as possible, so 
that the computational cost of the effective source does not prohibit its use in a numerical code. 

Fortunately, there are a two points which enable significant improvements: 

• Since the expansions are all about x, the coefficients of the polynomial do not change from grid point to grid 
point. They may change from one iteration to the next, however. 



* Even in 1 -|- ID and 2 -|- ID codes this is necessary because of the numerical integration involved. 
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• In some cases such as with circular orbits in Schwarzschild and Kerr spacetimes, this change between iterations 
is trivial and does not necessarily require recalculation of the effective source. 

This suggests an obvious optimization. The coefRcients are only computed once at the start of an iteration and then 
their numerical values are stored. The evaluation at each grid point then becomes simple multiplication by powers 
of (Ar, A0, A0); a relatively fast and computationally efficient operation. A further optimisation can be found by 
computing powers of (Ar, A0, A0) only once at the start of the simulation, providing the grid structure does not 
change. Altogether, this yields an enormous speed improvement - a factor of 50 — 100 in many cases. Similar tricks 
may also be employed with other parameters (mass, spin, etc.) which do not change through the lifetime of the 
simulation. Furthermore, if accuracy is important and delicate numerical cancellations are causing problems, this 
approach allows for the use of highly accurate methods such as Kahan jSCJj summation to minimize problems arising 
from numerical round-off. 

In addition to the numerical algorithm, it also important to consider the method for generating the code. Given 
the length of the expressions, it is impractical to manually type them in. Instead, we have directly generated C code 
from the Mathematica expressions and have made both available online [60j . 
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